
rm(list = ls(all = TRUE)) 

library(gdata)
library(reshape)
library(lubridate)
#library(erer)
library(readxl)   
library(stringr)
library(dplyr)
library(ggplot2)
library(gridExtra)
library(car)
library(foreign)
library(car)
library(miceadds)
library(jtools)
library(interplot)
#library(compactr)
library(lmtest)
library(sandwich)
library(multiwayvcov)
library(readxl)
library(data.table)
library(stringr)
library(haven)
library(ggrepel)
library(stargazer)
library(rstudioapi)


# set options
options(scipen=999)

# define a function to format names
simpleCap <- function(x) {
  s <- strsplit(x, " ")[[1]]
  paste(toupper(substring(s, 1,1)), substring(s, 2),
        sep="", collapse=" ")
}

# set working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#############################################################
## Function to calculate car for one country and one event
#############################################################

# requires df to have t and indices already calculated where t = 0 indicates event date
car.function <- function(est.win, ev.win, country, name, df, threshold, ff = FALSE){

# subset to country i in the estimation window
  est.sub <- df[df$t >= est.win[1] & df$t <= est.win[2] & df$Country.Name == country, ]

# remove na's
  est.sub <- est.sub[!is.na(est.sub$Spread) & !is.na(est.sub$Index_H) & !is.na(est.sub$Index_UM),]
  
# condition on at least one obs in estimation window
  if(nrow(est.sub) > 0){

# use ff indices if indicated 
  if(ff){
      fit <- lm(Spread ~ Index_H + Index_UM, data = est.sub) 
       } else{
    # regress spread on global index
    fit <- lm(Spread ~ Index, data = est.sub)  
}
  
# error message in case model produces non finite coefficients 
if(sum(is.na(fit$coefficients)) > 0){print(paste("Non-finite coefficients!", country, name, sep = " "))}
  
# condition on model fit 
if(is.finite(summary(fit)$adj.r.squared)){
 if(summary(fit)$adj.r.squared > threshold){
   
# subset to event window + estimation window for country i
    ev.sub <- df[((df$t >= est.win[1] & df$t <= est.win[2]) | (df$t >= ev.win[1] & df$t <= ev.win[2])) & 
                   df$Country.Name == country, ]

# predict spread for event and estimation windows given regression estimates
    if(ff == TRUE){
        ev.sub$Pred <- as.numeric(as.matrix(cbind(1, ev.sub$Index_H, ev.sub$Index_UM)) %*% fit$coefficients)
      }  else {
        ev.sub$Pred <- as.numeric(as.matrix(cbind(1, ev.sub$Index)) %*% fit$coefficients)
          }    
         
# calculate abnormal return
    ev.sub$ar <- ev.sub$Spread - ev.sub$Pred
    
# calculate standard deviation of abnormal return for estimation window only
    s.ar <- sqrt((1/(est.win[2] - est.win[1]))*sum((ev.sub$ar[ev.sub$t >= est.win[1] & ev.sub$t <= est.win[2]])^2))
    
# standard deviation of cumulative abnormal return 
    s.car <- sqrt(ev.win[2] - ev.win[1] + 1)*s.ar
    
# calculate cumulative abnormal return for event window
    car <- sum(ev.sub$ar[ev.sub$t >= ev.win[1] & ev.sub$t <= ev.win[2]])
    
# object to return
    car.sub <- c(as.character(country), # country name
                 car, s.car, round(car/s.car, digits = 3), # t stat
                 round(pt(-1*abs(car/s.car), ev.win[2] - ev.win[1]), digits = 3), # p-value
                 ifelse(pt(-1*abs(car/s.car), ev.win[2] - ev.win[1]) <= 0.05, '*', ''), # significance
                 summary(fit)$adj.r.squared)
    
    return(car.sub)
 } else {
   print(paste("Threshold failure!", country, name, sep = " "))
   return(NULL)
  }
  }
  } else{
  print(paste("No obs in estimation window!", country, name, sep = " "))
  return(NULL)
}
}

###############################################################################
## Function to calculate car for a set of countries given one event
###############################################################################

caar.function <- function(est.win, ev.win, event.date, region, name, data, threshold, ff = FALSE){

# create t variable where t = 0 indicates day of appointment  
  data$t <- as.integer(data$RDate - as.Date(as.character(event.date)))
  
# subset to event and estimation windows only
  df <- data[(data$t >= est.win[1] & data$t <= est.win[2]) | (data$t >= ev.win[1] & data$t <= ev.win[2]), ]

# identify treated units with data for both estimation and event windows
  treated1 <- unique(df$Country.Name[df$Region == region & df$t >= est.win[1] & df$t <= est.win[2]])
  treated2 <- unique(df$Country.Name[df$Region == region & df$t >= ev.win[1] & df$t <= ev.win[2]])
  countries <- treated1[treated1 %in% treated2]

# list control units available in estimation window
  for.index <- unique(df$Country.Name[!df$Country.Name %in% countries & df$t <= est.win[2]])
  
# calculate global index
  index.sub <- df[df$Country.Name %in% for.index, c('RDate', 'Country.Name', 'Spread')]
  index.sub <- index.sub %>%
    group_by(RDate) %>%
    summarize(Index = mean(Spread))
  df <- merge(df, index.sub, by = "RDate")
  
# calculate ff indices for each income group
  index.sub <- df[df$Country.Name %in% for.index, c('RDate', 'Country.Name', 'Spread', 'wb_class')]
  classes <- unique(index.sub$wb_class)
  for(g in 1:length(classes)){
    index.sub.g <- index.sub[index.sub$wb_class == classes[g], ]
    index.sub.g <- index.sub.g %>%
      group_by(RDate) %>%
      summarize('Index' = mean(Spread))
    names(index.sub.g)[2] <- paste('Index', classes[g], sep = '_')
    df <- merge(df, index.sub.g, by = "RDate", all.x = T) 
  }
  
# if any income groups missing create dummy vars so next line will run
if(!'Index_H' %in% names(df)){
  df$Index_H <- NA
}
if(!'Index_L' %in% names(df)){
    df$Index_L <- NA
}
  if(!'Index_UM' %in% names(df)){
    df$Index_UM <- NA
  }
  if(!'Index_LM' %in% names(df)){
    df$Index_LM <- NA
  }

# create income group matched index
df$Index_ff <- ifelse(df$wb_class == 'H', df$Index_H,
         ifelse(df$wb_class == 'L', df$Index_L,
                ifelse(df$wb_class == 'UM', df$Index_UM,
                       ifelse(df$wb_class == 'LM', df$Index_LM,
         NA))))

# drop any missing index values
df <- df[!is.na(df$Index_ff), ]

# create a holder to store results of car.function
  car.holder <- rep(NA, 7)

# condition on having data for at least one treated country available
if(length(countries) > 0){
  
# loop to calculate car for each country and add to holder
  for(i in 1:length(countries)){
    car <- car.function(est.win, ev.win, countries[i], name, df, threshold, ff)
    car.holder <- rbind(car.holder, car)
  } 
}
  
# error message in case no results
  if(length(na.omit(car.holder)) == 0){print(paste("No estimates for ", name, sep = ' '))}

# if no results return nulls
  if(length(na.omit(car.holder)) == 0){
    car.holder <- NULL
    caar.sub <- NULL
  } else {

# format car holder as data frame and format values
  car.holder <- as.data.frame(car.holder)
  car.holder <- na.omit(car.holder)
  car.holder[,2] <- as.numeric(as.character(car.holder[,2]))
  car.holder[,3] <- as.numeric(as.character(car.holder[,3]))

# calculate joint significance of cumulative abnormal returns
  caar <- mean(car.holder[,2])
  s.caar <- (1/nrow(car.holder)^2)*sqrt(sum((car.holder[,3])^2)) # see http://web.worldbank.org/archive/website01004/WEB/IMAGES/12852832.PDF
  
# create object to return  
  caar.sub <- c(as.character(name), caar, s.caar, 
                round(caar/s.caar, digits = 3), # t stat
                round(pt(-abs(caar/s.caar), ev.win[2] - ev.win[1]), digits = 2), # p value
                ifelse(pt(-abs(caar/s.caar), ev.win[2] - ev.win[1]) <= 0.05, '*', '')) # indicator of significance
  
# add staff name to cumulative.abnormal.returns
  car.holder$Staff.Name <- as.character(name)
  
}
  
  return(list(car.holder, caar.sub))
}



###########################################
## Function for multiple events
###########################################

es <- function(est.win, ev.win, events, data, threshold, ff = FALSE){
  # create holders
  car.holder <- rep(NA, 8)
  caar.holder <- rep(NA, 6)
  
  # loop across events
  for(t in 1:nrow(events)){
    event.holder <- NULL
    event.holder <- caar.function(est.win, ev.win, events[t,1], events[t,2], events[t,3], data, threshold, ff)
    if(length(event.holder[[1]]) > 0){
    car.holder <- rbind(car.holder, event.holder[[1]])
    caar.holder <- rbind(caar.holder, event.holder[[2]])
    }
  }
  
  # format results
  car.holder <- as.data.frame(car.holder)[-1,]
  names(car.holder) <- c("Country.Name", "car", "s.car", "t-stat", "p-value", "sig", "adj.r.sq", "Staff.Name")
  car.holder$car <- round(as.numeric(as.character(car.holder$car)), digits = 3)
  car.holder$s.car <- round(as.numeric(as.character(car.holder$s.car)), digits = 3)
  car.holder$adj.r.sq <- round(as.numeric(as.character(car.holder$adj.r.sq)), digits = 3)
  row.names(car.holder) <- NULL
  
  caar.holder <- as.data.frame(caar.holder)[-1,]
  names(caar.holder) <- c("Staff.Name", "caar", "s.caar", "t-stat", "p-value", "sig")
  caar.holder$Staff.Name <- as.character(caar.holder$Staff.Name)
  caar.holder$caar <- round(as.numeric(as.character(caar.holder$caar)), digits = 3)
  caar.holder$s.caar <- round(as.numeric(as.character(caar.holder$s.caar)), digits = 3)
  row.names(caar.holder) <- NULL
  
  wald <- sum(caar.holder$caar^2/caar.holder$s.caar^2) # note that this is a wald test not a t-test as used in the text
  p.wald <- 1-pchisq(wald, 1)
  sig.wald <- ifelse(1-pchisq(wald, 1) <= 0.05, '*', '')
  wald.test <- as.data.frame(t(c(round(wald, digits = 3), p.wald, sig.wald)))
  row.names(wald.test) <- NULL
  names(wald.test) <- c("wald.stat", "p.value", "sig")
 
  # return results
  AR <- list(car.holder, caar.holder, wald.test)
  names(AR) <- c('car', 'caar', 'wald.test')
  return(AR)
}





###########################################
## Function to Plot CAAR
###########################################

caar.plot <- function(es){
  # extract and count countries in each event
  car <- es[[1]]
  n.obs <- car %>% group_by(Staff.Name) %>% summarize(n.obs = n())
  
  # extract caar df and add n.obs
  coefs1 <- es[[2]]
  coefs1 <- merge(coefs1, n.obs, by = 'Staff.Name', all.x = T)
  
  # format caar df for plot
  coefs1$Staff.Name <- sapply(as.character(coefs1$Staff.Name), simpleCap)
  coefs1$caar <- as.numeric(as.character(coefs1$caar))
  coefs1$s.caar <- as.numeric(as.character(coefs1$s.caar))
  names(coefs1)[1] <- "Event"
  coefs1$Event <- factor(coefs1$Event, levels = unique(coefs1$Event[order(coefs1$caar)]))
  coefs1 <- coefs1[coefs1$n.obs > 1, ] 
  
  caar.plot <- ggplot(coefs1, aes(x = Event, y = caar)) +
    geom_hline(yintercept = 0, lty = 2, col = 'black') + geom_point(size = 1.5) +
    geom_errorbar(aes(x = Event, ymin = caar-qt(.975, n.obs)*s.caar, ymax = caar+qt(.975, n.obs)*s.caar), width=.5, size = .75) +
    labs(x = "", y = '') +
    theme_bw(base_size = 15) +  theme(legend.title = element_blank(), text = element_text(size=25)) + 
    coord_flip() 
  
  
  return(caar.plot)
  
  
}




###########################################
## Function to Make Quota Plot
###########################################

quota.plot <- function(es, qr){
  
  # format staff coefs
  car <- es[[1]]
  n.obs <- car %>% group_by(Staff.Name) %>% summarize(n.obs = n())
  
  coefs1 <- es[[2]]
  coefs1 <- merge(coefs1, n.obs, by = 'Staff.Name', all.x = T)
  coefs1 <- coefs1[coefs1$n.obs > 1, ]
  
  
  # format qr coefs
  qr.car <- qr[[1]]
  n.obs <- qr.car %>% group_by(Staff.Name) %>% summarize(n.obs = n())
  
  coefs2 <- qr[[2]]
  coefs2 <- merge(coefs2, n.obs, by = 'Staff.Name', all.x = T)
  coefs2 <- coefs2[coefs2$n.obs > 2, ]
  
  # misc formatting
  names(coefs1)[1] <- 'Staff.Name'
  coefs1$Staff.Name <- sapply(as.character(coefs1$Staff.Name), simpleCap)
  coefs1$caar <- as.numeric(as.character(coefs1$caar))
  coefs1$s.caar <- as.numeric(as.character(coefs1$s.caar))
  coefs2$caar <- as.numeric(as.character(coefs2$caar))
  coefs2$s.caar <- as.numeric(as.character(coefs2$s.caar))
  
  names(coefs2)[1] <- "Event"
  names(coefs1)[1] <- "Event"
  coefs1$Type <- "Staff Announcement"
  coefs2$Type <- "Quota Reform"
  
  coefs1 <- coefs1[order(coefs1$caar), ]
  coefs2 <- coefs2[order(coefs2$caar), ]
  
  coefs <- rbind(coefs2, coefs1)
  coefs$Event <- as.character(coefs$Event)
  
  coefs[1,1] <- 'IMF Reform 1'
  coefs[2,1] <- 'IMF Reform 2'
  coefs[3,1] <- 'DSK Arrested'
  
  coefs$Event <- factor(coefs$Event, coefs$Event)
  coefs$Type <- factor(coefs$Type, c("Staff Announcement", "Quota Reform"))
  
  
  
  qr.plot <- ggplot(coefs, aes(x = Event, y = caar, col = Type)) +
    geom_hline(yintercept = 0, lty = 2, col = 'black') + geom_point(size = 1.5) +
    scale_colour_manual(values = c('black', 'red')) +
    geom_errorbar(aes(x = Event, ymin = caar-qt(.975, n.obs)*s.caar, ymax = caar+qt(.975, n.obs)*s.caar), width=.5, size = .75) +
    labs(x = "", y = '') +
    theme_bw(base_size = 15) +  theme(legend.position = 'none', text = element_text(size=15)) + 
    coord_flip() 
  
  return(qr.plot)
  
}

###########################################
## Functions for Conditionality Analysis
###########################################

# scatterplot of caar vs conditionality
conditionality.scatter <- function(es, binding = T){
  e <- es[[2]]
  names(e)[1] <- 'Director'
  
  full <- merge(pred, e, by = 'Director')
  full$Correct1 <- ifelse(sign(full$Delta_BA1) == sign(full$caar), 'Yes', 'No')
  full$Correct2 <- ifelse(sign(full$Delta_BA2) == sign(full$caar), 'Yes', 'No')
  
  # scatter plot
  for(i in 1:nrow(full)){
    full$Director[i] <- simpleCap(full$Director[i])
  }
  full$Director <- ifelse(full$Director == 'Singh2', 'Singh (West)',
                          ifelse(full$Director == 'Singh', 'Singh (Asia)', full$Director))
  if(binding){
  plot <- ggplot(full, aes(x=caar, y=Delta_BA2)) +
    geom_point(size=3, shape=21, colour="lightgrey", aes(fill = factor(Correct2))) +
    scale_fill_manual(values=c("darkgrey", "red")) +
    geom_hline(yintercept = 0, linetype = 'dashed', color = 'grey') +
    geom_vline(xintercept = 0, linetype = 'dashed', color = 'grey') +
    geom_smooth(method=lm, se=T, linetype="dashed", color="darkgrey") +
    theme(legend.position="none") + 
    ylab('Change in Mean Binding Conditions') + xlab('Cumulative Average Abnormal Returns') +
    geom_label_repel(label = full$Director, segment.color = 'transparent', nudge_x = .2, nudge_y = 2)
  } else {
    plot <- ggplot(full, aes(x=caar, y=Delta_BA1)) +
      geom_point(size=3, shape=21, colour="lightgrey", aes(fill = factor(Correct2))) +
      scale_fill_manual(values=c("darkgrey", "red")) +
      geom_hline(yintercept = 0, linetype = 'dashed', color = 'grey') +
      geom_vline(xintercept = 0, linetype = 'dashed', color = 'grey') +
      geom_smooth(method=lm, se=T, linetype="dashed", color="darkgrey") +
      theme(legend.position="none") + 
      ylab('Change in Mean Conditions') + xlab('Cumulative Average Abnormal Returns') +
      geom_label_repel(label = full$Director, segment.color = 'transparent', nudge_x = .2, nudge_y = 2)
    
    
    
  }
  return(plot)
}
###########################################
## Function for Permutation Test
###########################################


perm.test <- function(ev, es, events, data, threshold, ff = T, reps){
  true.es <- es(ev, es, events, data, threshold, ff)
  holder <- as.data.frame(cbind(as.character(true.es[[2]][,1]), true.es[[2]][, 2]))
  names(holder) <- c('Staff.Name', 'True.CAAR')
  
  
  for(i in 1:reps){
    fake.events <- events
    fake.events$earliest.date <- fake.events$event.date - 60
    for(j in 1:nrow(fake.events)){
      fake.events$event.dates[j] <- sample(seq.Date(from = fake.events$earliest.date[j], to = fake.events$event.dates[j], 1), size = 1)
    }
    fake.es <- es(ev, es, fake.events, data, threshold, T)
    caar <- fake.es[[2]][,1:2]
    names(caar) <- c(names(caar)[1], paste('trial', i, sep = '.'))
    caar$Staff.Name <- as.character(caar$Staff.Name)
    holder <- merge(holder, caar, by = 'Staff.Name', all = T)
    fake.events <- NULL
  }
  
  holder[, 2] <- as.numeric(as.character(holder[, 2]))
  
  # results table              
  holder <- holder[!is.na(holder$True.CAAR),]
  holder$True.CAAR <- as.numeric(as.character(holder$True.CAAR))
  holder$p <- NA
  
  for(i in 1:nrow(holder)){
    p <- sum(abs(holder[i, 2]) < abs(holder[i, 3:(reps + 2)]), na.rm = T)/reps
    holder$p[i] <- p
    
    p <- NULL
    
  }
  
  # event (0,1)
  ptable <- holder
  ptable <- ptable[!is.na(ptable$True.CAAR), ]
  for(i in 1:nrow(ptable)){
    ptable$success[i] <- reps - sum(is.na(ptable[i,]), na.rm = T)
  }
  ptable <- ptable[,c('Staff.Name', 'p', 'success')]
  
  return(list(ptable, holder))
  
}

###########################################
## Data for Event Studies
###########################################

# first load and format wb classification data

class <- read_excel('_OGHIST.xlsx', sheet = "Country Analytical History")
names(class) <- class[5,]
names(class)[1:2] <- c('CountryCode', 'CountryName')
class <- class[c(11:239),]
class <- setDT(class)
long <- melt(class, id.vars = c("CountryCode","CountryName"), variable.name = "Year")

long$CountryName <- tolower(long$CountryName)
long$CountryName <- str_trim(long$CountryName, side = 'both')
long$CountryName <- ifelse(long$CountryName == 'bahamas, the', 'bahamas',
                           ifelse(long$CountryName == 'côte d\'ivoire', 'cote d\'ivoire',
                                  ifelse(long$CountryName == 'czechia', 'czech republic',
                                         ifelse(long$CountryName == 'egypt, arab rep.', 'egypt',
                                                ifelse(long$CountryName == 'hong kong sar, china', 'hong kong',
                                                       ifelse(long$CountryName == 'korea, rep.', 'korea',
                                                              ifelse(long$CountryName == 'lao pdr', 'lao people\'s dem. rep.',
                                                                     ifelse(long$CountryName == 'serbia and montenegro (former)', 'serbia and montenegro',
                                                                            ifelse(long$CountryName == 'türkiye', 'turkey',
                                                                                   ifelse(long$CountryName == 'united kingdom', 'uk',
                                                                                          ifelse(long$CountryName == 'venezuela, rb', 'venezuela',
                                                                                                 ifelse(long$CountryName == 'viet nam', 'vietnam', long$CountryName
                                                                                                 ))))))))))))
names(long)[4] <- 'wb_class'
names(long)[2] <- 'Country.Name'
long$Year <- as.integer(as.character(long$Year))

# load bonds data
data <- read.csv("_GFD_clean.csv")[,-1]
data$RDate <- as.Date(as.character(data$RDate))
# taking country europe out
dat <- data[data$Country.Name != "europe", ]
dat <- dat[dat$Country.Name != "world", ]

events <- read.csv('_events.csv')
events$event.dates <- as.Date(events$event.dates, format = '%m/%d/%y')

dat$Year <- year(dat$RDate)
long$Country.Name <- ifelse(long$Country.Name == 'taiwan, china', 'taiwan',
                            ifelse(long$Country.Name == 'tanzania', 'tanzania, united republic of',
                                   ifelse(long$Country.Name == 'uk', 'united kingdom',
                                          ifelse(long$Country.Name == 'vietnam', 'viet nam', long$Country.Name
                                                 ))))
dat <- merge(dat, long, by = c('Country.Name', 'Year'), all.x = T)

###########################################
## Data for Quota Reform
###########################################


treated.names <- tolower(c('Georgia', 'Honduras', 'Kenya', 'Kosovo', 'Serbia', 'Tunisia', 'Albania', 'Armenia',
                           'Cyprus', 'Greece', 'Jamaica', 'Pakistan', 'Seychelles', 'Ukraine', 'Colombia',
                           'Mexico', 'Poland', 'Morocco', 'Burkina Faso', 'Burundi', 'Chad', 'Ivory Coast',
                           'Ghana', 'Grenada', 'Guinea', 'Guinea-Bissau', 'Haiti', 'Kyrgyz', 'Liberia',
                           'Malawi', 'Mali', 'Niger', 'Sao Tome Principe', 'Sierra Leone', 'Solomon Islands',
                           'Yemen'))
dat.qr <- dat

dat.qr$Region <- ifelse(dat.qr$IMF.Program == 1, "imf", dat.qr$Region)

# set up list of events
event.dates <- c('2006-09-18', '2006-09-29', '2008-03-28', '2008-04-28', '2008-04-29', '2011-03-03',
                 '2015-12-16', '2011-05-14')
regions <- rep("imf", length(event.dates))
names <- event.dates
quota.events <- as.data.frame(cbind(event.dates, regions, names))


###########################################
## Data for Conditionality Analysis
###########################################

# load conditionality data
c <- read_excel('_IMFMonitor_Conditions_Raw.xlsx', sheet = 'Dataset')

# create indicators for binding and non-binding conditions then aggregate
c$BA1_indicator <- ifelse(c$`Condition Type` == 'IB' | 
                            c$`Condition Type` == 'QPC' |  
                            c$`Condition Type` == 'SB' | 
                            c$`Condition Type` == 'SPC' | 
                            c$`Condition Type` == 'PA' |
                            c$`Condition Type` == 'PC', 1, 0)

c$BA2_indicator <- ifelse(c$`Condition Type` == 'QPC' |  
                            c$`Condition Type` == 'SPC' | 
                            c$`Condition Type` == 'PA', 1, 0)


c.ag <- c %>%
  group_by(`Country Name`, `Country Code`, `Arrangement Date`, `Arrangement Type`, `Arrangement Amount`) %>%
  summarize(BA1 = sum(BA1_indicator), BA2 = sum(BA2_indicator), Waiver = sum(`Condition Waiver`, na.rm = T))

# add region variable
c.ag$`Country Name` <- tolower(c.ag$`Country Name`)

mcd <- c('afghanistan', 'algeria', 'armenia', 'azerbaijan', 'bahrain', 
         'djibouti', 'egypt', 'georgia', 'iran', 'iraq', 'jordan', 
         'kazakhstan', 'kuwait', 'kyrgyz republic', 'lebanon', 'libya', 'mauritania',
         'morocco', 'oman', 'pakistan', 'qatar', 'saudi arabia', 'somalia',
         'sudan', 'syria', 'tajikistan', 'tunisia', 'turkmenistan', 
         'united arab emirates', 'uzbekistan', 'yemen')

europe <- c('austria', 'belgium', 'cyprus', 'estonia', 'finland', 'france',
            'germany', 'greece', 'ireland', 'italy', 'latvia', 'lithuania',
            'luxembourg', 'malta', 'netherlands', 'portugal', 'slovak republic',
            'slovenia', 'spain', 'denmark', 'iceland', 'norway', 'sweden',
            'czechoslovakia', 'israel', 'switzerland', 'uk', 'hungary',
            'poland', 'bulgaria', 'croatia', 'romania', 'albania', 
            'bosnia and herzegovina', 'kosovo', 'macedonia', 
            'serbia and montenegro', 'belarus', 'moldova', 'ukraine', 
            'russian federation', 'turkey', 'north macedonia', 'yugoslavia', 'serbia',
            'czech republic')

# create region variable

af <- c('angola', 'benin', 'botswana', 'burkina faso', 'burundi', 
        'cabo verde', 'cameroon', 'central african republic', 'chad',
        'comoros', 'congo, dem. rep.', 'congo, rep.', 'cote d\'ivoire', 'equatorial guinea',
        'eritrea', 'ethiopia', 'gabon', 'gambia', 'ghana', 'guinea',
        'guinea-bissau', 'kenya', 'lesotho', 'liberia', 'madagascar',
        'malawi', 'mali', 'mauritius', 'mozambique', 'namibia',
        'niger', 'nigeria', 'rwanda', 'sao tome and principe', 'senegal',
        'seychelles', 'sierra leone', 'south africa', 'south sudan', 
        'swaziland', 'tanzania', 'togo', 'uganda', 'zambia', 'zimbabwe')


asia <- c('australia', 'new zealand', 'japan', 'hong kong', 'korea',
          'taiwan', 'singapore', 'bangladesh', 'brunei', 'cambodia',
          'china', 'india', 'indonesia', 'lao pdr',
          'malaysia', 'myanmar', 'mongolia', 'nepal', 'philippines',
          'sri lanka', 'thailand', 'vietnam', 'bhutan', 'fiji',
          'kiribati', 'maldives', 'marshall islands', 'micronesia',
          'nauru', 'palau', 'papua new guinea', 'samoa', 'solomon islands',
          'timor-leste', 'tonga', 'tuvalu', 'vanuatu')

west <- c('canada', 'mexico', 'us', 'argentina', 'bolivia', 'brazil',
          'chile', 'colombia', 'ecuador', 'guyana', 'paraguay',
          'peru', 'suriname', 'uruguay', 'venezuela', 'belize',
          'costa rica', 'el salvador', 'guatemala', 'honduras',
          'nicaragua', 'panama', 'antigua and barbuda', 'bahamas',
          'barbados', 'dominica', 'dominican republic', 'grenada',
          'haiti', 'jamaica', 'st. kitts and nevis', 'st. lucia',
          'st. vincent and the grenadines', 'trinidad and tobago',
          'cayman islands', 'bermuda')

c.ag$Region <- ifelse(c.ag$`Country Name` %in% mcd, 'mcd',
                        ifelse(c.ag$`Country Name` %in% af, 'af',
                               ifelse(c.ag$`Country Name` %in% west, 'west',
                                      ifelse(c.ag$`Country Name` %in% asia, 'asia', 
                                             ifelse(c.ag$`Country Name` %in% europe,'europe', NA)))))

# create director variable

events <- read.csv('_events.csv')
events$event.dates <- as.Date(events$event.dates, format = '%m/%d/%y')
events$start <- as.Date(events$start, format = '%m/%d/%y')
events$end <- as.Date(events$end, format = '%m/%d/%y')

c.ag$RDate <- as.Date(c.ag$`Arrangement Date`)
c.ag$Director <- NA

for(i in 1:nrow(c.ag)){
  sub <- events[events$start < c.ag$RDate[i] & 
                   events$end > c.ag$RDate[i] &
                   events$regions == c.ag$Region[i], ]
  if(nrow(sub) == 1){c.ag$Director[i] <- sub$names}
  sub <- NULL
}

# create year variable
c.ag$`Arrangement Amount` <- as.numeric(c.ag$`Arrangement Amount`)
c.ag$Year <- year(c.ag$`Arrangement Date`)

# summarize conditionality by director
d <- c.ag %>%
  group_by(Director, Region) %>%
  summarize(BA1 = mean(BA1), BA2 = mean(BA2), Waiver = mean(Waiver), Loan = mean(`Arrangement Amount`, na.rm = T))

d <- d[d$Director != 'gondwe', ]

names(events)[3] <- 'Director'
d <- merge(events, d, by = 'Director')
names(events)[3] <- 'names'


# add predecessor variables

pred <- d

pred$Pred <- ifelse(pred$Director == 'horiguchi', 'neiss',
                    ifelse(pred$Director == 'burton', 'horiguchi',
                           ifelse(pred$Director == 'singh', 'burton',
                                  ifelse(pred$Director == 'rhee', 'singh',
                                         ifelse(pred$Director == 'belka', 'deppler',
                                                ifelse(pred$Director == 'borges', 'belka',
                                                       ifelse(pred$Director == 'moghadam', 'borges',
                                                              ifelse(pred$Director == 'thomsen', 'moghadam',
                                                                     ifelse(pred$Director == 'khan', 'abed',
                                                                            ifelse(pred$Director == 'azour', 'ahmed',
                                                                                   ifelse(pred$Director == 'singh2', 'loser',
                                                                                          ifelse(pred$Director == 'eyzaguirre', 'singh2',
                                                                                                 ifelse(pred$Director == 'werner', 'eyzaguirre', NA)))))))))))))

names(d) <- c('Pred', 'Event.Date', 'regions', 'Start', 'End', 'Region', 'Pred_BA1', 'Pred_BA2', 'Pred_Waiver', 'Pred_Loan')

pred <- merge(pred, d[,c("Pred", "Pred_BA1", "Pred_BA2")], by = 'Pred')

# calculate deltas

pred$Delta_BA1 <- pred$BA1 - pred$Pred_BA1
pred$Delta_BA2 <- pred$BA2 - pred$Pred_BA2

###########################################
## Main Results
###########################################

# event study models
es1 <- es(c(-180, -5), c(0, 1), events, dat, .25, T)
es2 <- es(c(-180, -5), c(0, 3), events, dat, .25, T)
es3 <- es(c(-180, -5), c(0, 5), events, dat, .25, T)

# caar plots
fig_1a <- caar.plot(es1)
fig_1b <- caar.plot(es2)
fig_1c <- caar.plot(es3)

ggsave(filename = 'Figures/fig_1a.png', plot = fig_1a,
       width = 1500, height = 1500, units = "px")

ggsave(filename = 'Figures/fig_1b.png', plot = fig_1b,
       width = 1500, height = 1500, units = "px")

ggsave(filename = 'Figures/fig_1c.png', plot = fig_1c,
       width = 1500, height = 1500, units = "px")


# quota reform analysis and plot
qr <- es(c(-180, -5), c(0, 3), quota.events, dat.qr, 0, T)
fig_2a <- quota.plot(es2, qr)

ggsave(filename = 'Figures/fig_2a.png', plot = fig_2a,
       width = 2000, height = 1500, units = "px")


# conditionality plot
fig_2b <- conditionality.scatter(es3, F)
fig_a5 <- conditionality.scatter(es3, T)

ggsave(filename = 'Figures/fig_2b.png', plot = fig_2b,
       width = 1800, height = 1500, units = "px")
ggsave(filename = 'Figures/fig_a5.png', plot = fig_a5,
       width = 2200, height = 1500, units = "px")

###########################################
## Identity Characteristics + Principal Heterogeneity
###########################################

# load unga voting data and subset to useful vars
un <- read.csv('_IdealpointestimatesAll_Jun2024.csv')
un <- un[, c('IdealPointAll', 'Countryname', 'session')]

# create year variable and g-5 dummy
un$year <- un$session + 1945
un$g5 <- ifelse(un$Countryname %in% c('United States', 'United Kingdom',
                                      'German Federal Republic', 'France',
                                      'Japan'), 1, 0)
# subset to g-5
un <- un[un$g5 == 1, ]

# create pairwise dataframe and rename
un.dummy <- un
names(un) <- c('IdealPoint_A', 'Countryname_A', 'session', 'year', 'g5')
names(un.dummy) <- c('IdealPoint_B', 'Countryname_B', 'session', 'year', 'g5')

un <- merge(un, un.dummy, by = c('year', 'session', 'g5'),)
un <- un[un$Countryname_A != un$Countryname_B, ]

# create distance measure
un$IdealPoint_Dif <- abs(un$IdealPoint_A - un$IdealPoint_B)

# average by year
un <- un %>%
  group_by(year) %>%
  summarize(PHet = mean(IdealPoint_Dif))

# extract event study estimates and merge in years, prinicpal heterogeneity
e <- es1[[2]]
y <- events[, c('event.dates', 'names')]
y$year <- year(y$event.dates)
names(y)[2] <- c('Staff.Name')

e <- merge(e, y, by = 'Staff.Name', all.x = T)
e <- merge(e, un, by = 'year')

# format staff names
for(i in 1:nrow(e)){
e$Staff.Name[i] <- simpleCap(e$Staff.Name[i])
}

e$Staff.Name <- ifelse(e$Staff.Name == 'Singh2', 'Singh (West)',
                       ifelse(e$Staff.Name == 'Singh', 'Singh (Asia)', e$Staff.Name))
       
# dummy for global south                                                                                        ifelse(e$Staff.Name == '', '',)
e$gs <- ifelse(e$Staff.Name %in% c('Thomsen', 'Burton', 'Belka', 'Borges'), 0, 1)


# dummy for neoliberal phd
e$phd <- ifelse(e$Staff.Name %in% c('Loser', 'Horiguchi', 'Khan',
                                    'Singh (West)', 'Singh (Asia)',
                                    'Belka', 'Rhee', 'Azour',
                                    'Burton', 'Werner', 'Moghadam', 'Borges'), 1, 0)

# histogram by education
fig_3a <- ggplot(e, aes(x=caar, fill=factor(phd))) + 
  scale_fill_manual(labels = c("No", "Yes"), values = c('red', 'orange')) +
  geom_density(alpha=.2) + 
  geom_vline(xintercept = mean(e$caar[e$phd == 1]), linetype = 'dashed', color = 'black') +
  geom_vline(xintercept = mean(e$caar[e$phd == 0]), linetype = 'dashed', color = 'black') +
  labs(x="Cumulative Average Abnormal Returns", y="Density", fill = 'Economics Ph.D.')

# histogram by gs
fig_3b <- ggplot(e, aes(x=caar, fill=factor(gs))) + 
  scale_fill_manual(labels = c("No", "Yes"), values = c('purple', 'blue')) +
  geom_density(alpha=.2) + 
  geom_vline(xintercept = mean(e$caar[e$gs == 1]), linetype = 'dashed', color = 'black') +
  geom_vline(xintercept = mean(e$caar[e$gs == 0]), linetype = 'dashed', color = 'black') +
  labs(x="Cumulative Average Abnormal Returns", y="Density", fill = 'Global South')
  
# plot by phet
fig_a6 <- ggplot(e, aes(x=caar, y=exp(PHet))) +
  geom_point(size=3, shape=23, colour="black", fill = 'black') +
  geom_vline(xintercept = 0, linetype = 'dashed', color = 'grey') +
  geom_smooth(method=lm, se=T, linetype="dashed", color="darkgrey") +
  theme(legend.position="none") + 
  ylab('Principal Heterogeneity') + xlab('Cumulative Average Abnormal Returns') +
  geom_label_repel(label = e$Staff.Name, segment.color = 'transparent', nudge_x = .15, nudge_y = .025)

ggsave(filename = 'Figures/fig_3a.png', plot = fig_3a,
       width = 2000, height = 1500, units = "px")
ggsave(filename = 'Figures/fig_3b.png', plot = fig_3b,
       width = 2000, height = 1500, units = "px")
ggsave(filename = 'Figures/fig_a6.png', plot = fig_a6,
       width = 2200, height = 1500, units = "px")



# extract country level results and format names
c <- es1[[1]]
for(i in 1:nrow(c)){
  c$Staff.Name[i] <- simpleCap(c$Staff.Name[i])
  
}

c$Staff.Name <- ifelse(c$Staff.Name == 'Singh', 'Singh (Asia)',
                       ifelse(c$Staff.Name == 'Singh2', 'Singh (West)', c$Staff.Name))

# merge in characteristic indicators
c <- merge(c, e[, c('Staff.Name', 'phd', 'gs')], by = 'Staff.Name', all.x = T)

# regress abnormal returns on identity indicators
fit <- summary(lm(car ~ phd + gs, data = c))

stargazer(fit$coefficients, summary = FALSE, 
          star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c("Economics Ph.D.", "National of Global South"),
          dep.var.labels = "Cumulative Abnormal Returns", 
          out="./Tables/table_a3.tex",
          label="table_a3")


###########################################
## Plot Conditionality Over Time
###########################################

t <- events[, 2:5]
names(t) <- c('Region', 'Director', 'Start', 'End')
t$Start <- as.Date(t$Start, format = "%m/%d/%y")
t$End <- as.Date(t$End, format = "%m/%d/%y")

t$Start_Year <- year(t$Start)
t$End_Year <- year(t$End)
t <- t[, c(1,2,5,6)]
t <- merge(t, pred[, c("Director", "BA1", "BA2")], by = "Director")


df <- rep(NA, 6)
for(i in 1:nrow(t)){
  l <- length(t$Start_Year[i]:t$End_Year[i]) - 1
  
  h <- cbind(rep(t$Director[i], l), rep(t$Region[i], l), rep(t$BA1[i], l),
             rep(t$BA2[i], l), (t$Start_Year[i] + 1):t$End_Year[i], c(simpleCap(t$Director[i]), rep(NA, l-1)))
  
  df <- rbind(df, h)  
  
}

df <- as.data.frame(df)
names(df) <- c("Director", "Region", "BA1", "BA2", "Year", "Label")
df$BA1 <- as.numeric(as.character(df$BA1))
df$BA2 <- as.numeric(as.character(df$BA2))

df <- df[!is.na(df$Region),]

df$Label <- ifelse(df$Label == 'Singh2', 'Singh', df$Label)
df$Year <- as.numeric(as.character(df$Year))
for(i in 1:nrow(df)){
  df$Region[i] <- simpleCap(df$Region[i])
}

fig_a4a <- ggplot(aes(x=Year, y=BA1, colour = Director, group = Region), data = df) + 
  geom_line() + geom_point() + ylab('Mean Conditions Imposed') +
  geom_label_repel(label=df$Label, nudge_x = .5, nudge_y = 3, segment.color = 'transparent') + 
  theme(legend.position="none", axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  xlim(2000, 2024) +
  facet_wrap(~Region) 

fig_a4b <- ggplot(aes(x=Year, y=BA2, colour = Director, group = Region), data = df) + 
  geom_line() + geom_point() + ylab('Mean (Binding) Conditions Imposed') +
  geom_label_repel(label=df$Label, nudge_x = .5, nudge_y = 3, segment.color = 'transparent') + 
  theme(legend.position="none", axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  xlim(2000, 2024) +
  facet_wrap(~Region) 

ggsave(filename = 'Figures/fig_a4a.png', plot = fig_a4a,
       width = 3000, height = 1500, units = "px")
ggsave(filename = 'Figures/fig_a4b.png', plot = fig_a4b,
       width = 3000, height = 1500, units = "px")



###########################################
## Permutation Test and Table
###########################################

set.seed(2000)
# run permutation tests
ptest1 <- perm.test(c(-180, -5), c(0,1), events, dat, .25, T, 1000)
ptest2 <- perm.test(c(-180, -5), c(0,5), events, dat, .25, T, 1000)

# extract test results
ptable1 <- ptest1[[1]]
ptable2 <- ptest2[[1]]

# merge test results
ptable <- merge(ptable1, ptable2, by = c("Staff.Name"), all.x = TRUE, all.y = TRUE)

ptable$Staff.Name.Upper <- toupper(ptable$Staff.Name)
ptable$p.x <- ptable$p.x
ptable$p.y <- ptable$p.y

ptable <- ptable[, c("Staff.Name.Upper", "p.x", "p.y")]
# create table
stargazer(ptable[order(ptable[, 2], decreasing = FALSE),], summary = FALSE, 
          keep = c("Staff.Name.Upper", "p.x", "p.y"),
          out = "./Tables/table_a2.tex", 
          label = "table_a2")


